Competing states in the SU(3) Heisenberg model on the honeycomb lattice: 
Plaquette valence-bond crystal versus dimerized color-ordered state 
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Conflicting predictions have been made for the ground state of the SU(3) Heisenberg model on the honeycomb 
lattice: Tensor network simulations found a plaquette order [Zhao et al, Phys. Rev. B 85, 134416 (2012)], where 
singlets are formed on hexagons, while linear flavor- wave theory (LFWT) suggested a dimerized, color ordered 
state [Lee and Yang, Phys. Rev. B 85, 100402 (2012)]. In this work we show that the former state is the true 
ground state by a systematic study with infinite projected-entangled pair states (iPEPS), for which the accuracy 
can be systematically controlled by the so-called bond dimension D. Both competing states can be reproduced 
with iPEPS by using different unit cell sizes. For small D the dimer state has a lower variational energy than 
the plaquette state, however, for large D it is the latter which becomes energetically favorable. The plaquette 
formation is also confirmed by exact diagonalizations and variational Monte Carlo studies, according to which 
both the dimerized and plaquette states are non-chiral flux states. 



PACS numbers: 67.85.-d, 71.10.Fd, 75.10.Jm, 02.70.-c 

I. INTRODUCTION 

Systems of fermions with multiple flavors have attracted 
increasing interest recently thanks to the proposals to experi- 
mentally realize such systems with ultra-cold fermionic atoms 
in optical lattices ^""^ and the rapid experimental progress in the 
field.^"^^ In general these systems can be described by a Hub- 
bard model with N flavors (or colors) of fermions. In the limit 
of strong on-site repulsion and an integer filling per lattice site, 
the system is in a IMott insulating state, and the low-energy 
physics is captured by the SU(N) Heisenberg model. These 
models give rise to a rich variety of exotic quantum states, 
such as different Neel-type states,^^"^^ generalized valence- 
bond solids,^^"^^ algebraic N-flavor liquids,^^"^^ chiral spin 
liquids,^^'^^'^^'^"^, and more.^^~^^ 

The SU(N) Heisenberg Hamiltonian is given by 

where the first sum goes over nearest neighbors pairs and a, 
13 run over the N possible colors (flavors) at each site. Pij is a 
permutation operator which exchanges colors on neighboring 
sites. In the present work we will focus on the case of one 
particle per site (corresponding to the fundamental represen- 
tation) and TV = 3 on the honeycomb lattice. 

In general the theoretical study of these models is very chal- 
lenging, particularly because in many cases these models ex- 
hibit a negative sign problem in Quantum IMonte Carlo sim- 
ulations, in contrast to the A/" = 2 case on bipartite lattices. 
Therefore, other methods have to be used to study the physics 
of these models. ]V[ean-field theory typically fails to correctly 
predict the ground state. In most cases the classical solution 
exhibits an infinite degeneracy, which is lifted upon including 
quantum fluctuations. 



One way to go beyond the simple mean-field (Hartree) solu- 
tion is through linear flavor- wave theory, which takes into ac- 
count quantum fluctuations on top of a Hartree solution at the 
harmonic level. This method has successfully accounted for 
the three- sublattice order in the SU(3) Heisenberg model on 
the square lattice. Further quantum fluctuations can in prin- 
ciple be included by taking higher-order terms into account, 
but this has not been achieved yet. 

A powerful class of methods which enable a systematic 
study of the solution upon adding quantum fluctuations are 
tensor network algorithms. The most famous method in this 
class is the density matrix renormalization group (DIMRG) 
method,^^ which is the state-of-the-art method to simulate 
(quasi-) one dimensional systems. DIMRG is based on a vari- 
ational ansatz called matrix product state (]V[PS), where the 
coefficients of the wave function are efficiently encoded by 
a product of matrices. Substantial progress has also been 
made in the simulation of two-dimensional systems with ex- 
tensions of the IMPS to higher dimensions, a so-called pro- 
jected entangled-pair state (PEPS) or tensor product state. ^^'"^^ 
As in an IMPS, the accuracy of a PEPS can be controlled by 
the so-called bond dimension D. 

Based on these two approaches, conflicting results have 
been reported for the SU(3) Heisenberg model on the hon- 
eycomb lattice recently."^ ^'"^^ In both studies the bilinear- 
biquadratic spin-1 model has been considered, which includes 
the SU(3) Heisenberg model as a special point in the phase 
diagram (when the coefficients of the bilinear and biquadratic 
terms are equal and positive). LEWT"^^ predicts a dimerized, 
color-ordered state in a 18-site unit cell, depicted in Fig. 1, 
whereas in the variational tensor network study in Ref. 42 
a plaquette order has been found by using a 6-site unit cell. 
LFWT is not a variational method and therefore the energies 
from the two approaches cannot be directly compared. Fur- 
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FIG. 1: (Color online) Schematic illustration of the two compet- 
ing states of the SU(3) Heisenberg model on the honeycomb lat- 
tice. Thick (thin) bonds correspond to low (high) bond energies, 
(a) Color-ordered, dimerized state obtained with a 18-site unit cell 
(shaded in grey), (b) Plaquette state with a 6- site unit cell which has 
no color order. 



thermore, a 18-site unit cell, which would be compatible with 
the dimerized state, has not been tested in the tensor network 
study, and thus it is still an open question what the true ground 
state is. 

The aim of our work is to unambiguously identify the true 
ground state of the SU(3) Heisenberg model on the honey- 
comb lattice by a systematic study of the energetics of the two 
competing states by means of infinite PEPS (iPEPS) with dif- 
ferent unit cell sizes and bond dimensions. We show that the 
dimerized state predicted by LFWT can be understood as a 
low-entanglement solution which is reproduced with iPEPS 
for small bond dimensions. This state, however, is metastable 
for large bond dimensions, and the true ground state is the pla- 
quette state found in the previous tensor network study. We 
provide further support for this state also from exact diago- 
nalization results up to system size Ng = 24. Finally, us- 
ing Gutz wilier-projected free-fermion wavefunction, we char- 
acterize the competing states based on the properties of the 



fermionic wave function. 

The paper is organized as follows: In Sec. II A we provide a 
brief introduction to iPEPS, in Sec. II B we present the iPEPS 
results, in Sec. Ill the ED results and in Sec. IV the varia- 
tional Monte Carlo (VMC) results from Gutzwiller projected 
fermionic wavef unctions. Finally we summarize our conclu- 
sions in Sec. V. 



II. INFINITE PROJECTED ENTANGLED-PAIR STATE 
(IPEPS) 

A. Method 

In this section we give a short summary of the tensor net- 
work method used in this work, and point out the differences 
to the approach used in Ref. 42. For further details on the 
method we refer to previous works,"^^'"^^'"^"^ in particular also 
Ref. 31 where we used the same approach for the SU(4) 
Heisenberg model on the honeycomb lattice. 

Each tensor network algorithm has three essential ingredi- 
ents: (1) the structure of the tensor network ansatz, (2) the op- 
timization method (i.e. how to find the optimal values for the 
tensors to have an approximate representation of the ground 
state), and (3) the method used to compute expectation values 
of observables (i.e. the contraction of the tensor network). 

(1) The tensor network ansatz we use is a projected 
entangled-pair state (PEPS).^^'"^^ It is a variational ansatz 
aimed at efficiently representing ground states of two- 
dimensional lattice models. The coefficients of the wave func- 
tion are obtained by taking the trace of a product of tensors, 
with one tensor per lattice site. Each tensor has a physical in- 
dex which carries the local Hilbert space of a lattice site of di- 
mension d, and auxiliary indices which connect to the neigh- 
boring tensors on the lattice. These auxiliary indices have a 
certain dimension D which is called the bond dimension. On 
the square lattice each tensor T^^^^ has dD"^ elements, whereas 
on the honeycomb lattice each tensor T^^^ has dD^ elements. 
The numbers stored in these tensors are the variational pa- 
rameters of the ansatz, i.e. the larger D the more variational 
parameters, and therefore the (potentially) more accurate the 
ansatz. The special case of D = 1 corresponds to a prod- 
uct state with a vector on each site. An infinite PEPS 
(iPEPS)^^ consists of a unit cell of different tensors, which is 
periodically repeated on the lattice to represent a state directly 
in the thermodynamic limit. Using different unit cell sizes en- 
ables to represent states with different types of translational 
symmetry breaking. 

(2) An approximate representation of the ground state is 
found by performing an imaginary time evolution of a ran- 
dom initial iPEPS. The imaginary time evolution operator 
ex.-p{—/3H) is split into a product of two-body operators via 
a second order Trotter-Suzuki decomposition (see Ref. 44). 
Multiplying such a two-body operator to the iPEPS increases 
the bond dimension of the corresponding bond between the 
sites the operator is acting on. To constrain the computational 
cost the auxiliary space of the bond has to be truncated down 
to the original bond dimension D after a two-body operator 
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FIG. 2: (Color online) Local order parameters obtained with iPEPS 
obtained with two different unit cell sizes. The thickness of the bonds 
is proportional to the square of the energy of the corresponding bond. 
Each pie chart on a site visualizes the local color density, (a) The 
color-ordered, dimerized state obtained with a 18-site unit cell for 
D = 2. (b) Plaquette state which has no color order: Each color has 
the same density on each site (D = 16). 



has been applied. There are different ways to perform this 
truncation, as discussed e.g. in Ref 44. With the full update a 
bond is truncated in an optimal way by taking into account the 
whole wave function to find the relevant subspace. The simple 
update^^'^^ is computationally cheaper since it involves only 
local tensors surrounding the bond to be truncated, however, 
it is not optimal. In the present work we use the more accurate 
full update (in contrast to the previous study in Ref. 42). 

(3) To evaluate observables the tensor network needs to 
be contracted (by computing the trace of the product of all 
tensors). This contraction can only be done approximately 
in polynomial time. As in previous works we use the cor- 
ner transfer matrix method"^^'"^^ generalized to arbitrary unit 
cell sizes. We map the honeycomb lattice onto a brick- wall 
square lattice as explained in Ref. 3 1 . We checked that quan- 
tities of interest are converged in the "boundary" dimension 
X, which controls the accuracy of the contraction. 

In order to reduce the computational cost we use tensors 
with Ijq symmetry, which is a discrete subgroup of SU(3).^^'^^ 
The tensors then acquire a block structure, similar to a block 
diagonal matrix. 




1/D or 1/N 



FIG. 3: (Color online) iPEPS simulation results for the variational 
energy as a function of inverse bond dimension 1/D for different unit 
cell sizes, compared to the exact diagonalization (ED) results and 
the lowest energy states from variational Monte Carlo (VMC). For 
bond dimensions D < 8 the color-ordered, dimerized state has the 
lowest variational energy, for larger D the plaquette state becomes 
energetically favorable. The dotted lines are a guide to the eye. 



B. IPEPS results 

Here we present a systematic study of the solution for the 
ground state as a function of the bond dimension D in iPEPS, 
which controls the accuracy of the ansatz, and also the amount 
of quantum fluctuations (or entanglement) taken into account, 
as we explain in the following. We consider results for differ- 
ent unit cell sizes: 2x2, and the 6-site and the 18-site unit 
cell shaded in grey in Fig. l(a-b). 

A D = 1 iPEPS corresponds to a product state (a site- 
factorized wave function), i.e. a non-entangled state. The 
energy per site is = 0, which can be easily verified an- 
alytically. The state is infinitely degenerate: for example all 
possible coverings where two nearest-neighbor sites exhibit 
different colors have the same energy (or more generally, the 
energy is minimized if the states on neighboring sites are or- 
thogonal). 

With D = 2 short-range quantum fluctuations are taken 
into account. It turns out that they lift the infinite degeneracy. 
If a 18-site unit cell is used in iPEPS, the quantum fluctua- 
tions select the same state as predicted by LFWT"^^ : the dimer- 
ized, color-ordered state shown in Fig. 1(a). In Fig. 2(a) we 
visualize different local quantities obtained from the iPEPS: 
The thickness of the bonds is proportional to the square of the 
corresponding bond energy, and the pie charts show the local 
color density of each color. On each site one of the colors 
is dominant, and the pattern of the dominant colors matches 
the one shown in Fig. 1(a). The state is clearly dimerized: 
the two sites in each dimer have different colors (e.g. green 
and red), and each dimer is surrounded by four sites where 
the third color is dominant (e.g. blue). The energy per site is 
Ed=2 = -0.553. 

If we take smaller unit cells this state cannot be represented, 
and the variational energy is higher, e.g. —0.421 in the 2 x 2 
unit cell and —0.362 in the 6-site unit cell (using Zq symmet- 
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ric tensors). 

By further increasing D the variational energies in the dif- 
ferent unit cells are lowered, as shown in Fig. 3. The energy 
obtained with the 6- site unit cell decreases more rapidly with 
D than the energy obtained with the 18-site unit cell, and it 
becomes lower for D > 10. The state obtained in the 6- site 
unit cell is the plaquette state, shown in Figs. 1(b) and 2(b), 
where low energy bonds are formed around hexagons. On all 
sites each color has the same density, i.e. the state does not ex- 
hibit color-order. This state has already been found in Ref. 42. 
So, the dimerized state is not the true ground state, but only 
a metastable state which appears when some, but not all of 
the quantum fluctuations are taken into account. We can call 
it a low-entanglement solution, since it is energetically favor- 
able for small values of D (e.g. D = 2), at which the iPEPS 
represents only a slightly entangled state. 

By extrapolating the energy of the plaquette state to the in- 
finite D limit, we expect the ground state energy to lie in be- 
tween —0.723 and —0.720. From the slopes of the curves we 
do not expect another crossing at larger D. 

In Fig. 4(a) we show results for the local ordered moment 
m, given by 




where = |q^)(/3| are the generators of SU(3) and a, /3 run 
over all flavor indices. For the plaquette state m vanishes for 
D > 10 (i.e. there is no color order). For the dimerized state 
m remains finite for all values of D, however, m is strongly 
suppressed with increasing D. Extrapolating minl/D yields 
a finite value, however, since the extrapolated value is small it 
is difficult to conclude whether the SU(3) symmetry is broken 
in the infinite D limit or not. 

In Fig. 4(b) we plot the difference between the highest bond 
energy and the lowest bond energy in the unit cell, 

AE = max(£;5) - min(^5), (3) 

which measures the strength of the dimerization or plaquette 
order. The plaquette order is suppressed with increasing D, 
however, it seems to tend to a finite value in the infinite D 
limit, which shows that the ground state indeed has long-range 
plaquette order in the thermodynamic limit. 

We conclude this section with a remark: The plaquette state 
is compatible with the 18-site unit cell, thus the simulations 
with this unit cell for I) > 8 should in principle yield the 
plaquette state. The reason why we remain in the dimer state 
for I) > 8 is due to metastability. Since the states for lower 
D are used as an initial state for simulations at larger D one 
has to overcome an energy barrier to get from the dimer state 
into the plaquette state when moving from I^ = 8tol^ = 10. 
This does not seem to occur, at least not in the simulated time 
scales, i.e. the simulation is stuck in the metastable dimer 
state. We can exploit this fact to compare the energies of the 
two states at large D. 




0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 
1/D 1/D 



FIG. 4: (Color online) iPEPS simulation results of order parame- 
ters as a function of inverse bond dimension 1/D different unit cell 
sizes. Dotted lines are a guide to the eye. (a) Local moment which is 
suppressed with increasing bond dimension. For the plaquette state 
it vanishes for D > 10, indicating absence of color-order. Lower 
right panel: Difference in energy between the strongest bond and 
the weakest bond in the unit cell, showing that the plaquette state is 
stable in the infinite D limit. 



III. ED RESULTS 

In order to corroborate the iPEPS findings we also per- 
formed exact diagonalizations of the SU(3) Heisenberg model 
on two clusters consisting of Ns = 18 and Ns = 24 sites. 
Both clusters are compatible with the plaquette state, how- 
ever the Ns = 18 is particular in that it contains additional 
loops of length six which wrap around the torus. As a con- 
sequence this cluster allows for more "plaquette" states than 
clusters with larger circumference. Regarding the magnetic 
dimer state, only the Ng = 18 cluster is compatible with this 
type of order. 

We first calculated the ground state energy per site, which is 
shown together with the energies of the iPEPS and the VMC 
approaches in Fig. 3. The energy per site of the Ng = IS clus- 
ter is -0.76928, while the value for Ns = 24 is -0.72687, 
very close to the large D value inferred from the iPEPS re- 
sults. 

Then, we determined the low energy spectrum resolved ac- 
cording to spatial quantum numbers and for different mag- 
netization sectors, i.e. different values of the three colors 
{Na^ Nbj Nc)- These spectra are displayed in Fig. 5. The 
spectrum of the A^^ = 18 cluster (top panel) shows a signifi- 
cant gap between the non-degenerate ground state and a set of 
states at an energy of about ~ —12.6 (encircled levels). The 
two levels highlighted with full line circles are precisely the 
two levels which are expected for a non-magnetic plaquette 
state with a three-fold ground state degeneracy. The (two-fold 
degenerate) level highlighted with a dashed circle is an artifact 
of the enhanced spatial symmetry of the Ns = IS honeycomb 
cluster. Finally the magnetic level encircled with a dotted cir- 
cle would belong to a hypothetical tower of state structure for 
the magnetic ordering of the dimer state. However since there 
is no clear ordering between the levels making up the plaque- 
tte state or the dimer state, the A^^ = 18 cluster does not seem 
to be helpful in deciding between the two competing states. 

The situation is much more clear on the A^g = 24 cluster. 
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(Color online) Real space bond energy correlation 
- {Pij){Pki) obtained from ED on clusters of (a) Ns — 18 
and (b) Ns = 24 sites. The value of the correlation is proportional to 
the plotted bond strength. Positive (negative) correlations are shown 
as full blue (red, dashed) lines. The reference bond is denoted as the 
bold black bond. 
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FIG. 5: (Color online) Low energy spectra from ED for Ns — 1^ (top 
figure) and Ns = 24 (bottom figure). See main text for explanation 
of the different symbols. 



where there is an almost perfect degeneracy between the non- 
degenerate ground state (F Al representation) and the twofold 
degenerate level in the sector i^T Al, fully consistent with a 
plaquette-like spatial symmetry breaking. 

We have also calculated the color-color correlation func- 
tions in the ground states of both clusters (not shown), and 
even though the correlations on the A^s = 18 cluster are con- 
sistent in sign with the dimerized magnetic state, the actual 
correlations are quite weak beyond the nearest neighbor corre- 
lations, and thus do not lend strong support for a magnetically 
ordered state. 

Finally we show the real-space bond-energy correlations in 
Fig. 6. Common to both clusters is the fact that these correla- 
tions extend through the entire cluster, i.e. the correlations are 
quite long ranged. When it comes to the signs of the corre- 
lations, it turns out that the A^g = 18 correlations suffer from 
an admixture of the additional plaquette states, which are an 
artifact of this cluster. The correlations in the ground state of 
the Ns = 24 cluster do not suffer from this and most bonds 
match the expectations of a plaquette ordered singlet state (see 
e.g. Ref. 52 for a discussion of the real-space correlations of 
SU{2) plaquette states on the honeycomb lattice). 

In conclusion the ED results are consistent with the pla- 
quette ordering scenario put forward by the Tensor Network 
and the iPEPS approach, even though they cannot strictly rule 



out the magnetic ordering scenario due to the lack of larger 
clusters which are compatible with both competing states, and 
which would allow an unbiased comparison. 



IV. VARIATIONAL MONTE CARLO 

A. Introduction 

The fermionic representation of the different colors has 
proven to be a valuable tool to understand the ground state 
properties of the SU(2) and SU(N) Heisenberg models. The 
permutation operator Pij that exchanges the colors between 
sites i and j can be written as 

N 

PiJ = - E /I«40-)//3«/a(i) (4) 

a,/3=l 

in the fermionic representation, where f^{i) creates a fermion 
of color a at site i, and fd^) annihilates it. To treat the four 
fermion term, it is customary to introduce a bond mean-field 
decoupling of the Pij of the form 

(5) 

so that = ^(^i j) Pf^J describes free fermions. The hop- 
ping amplitude Xi,j is determined self-consistently as the ex- 
pectation value 



N 



(6) 



OL=l 



taken in the ground state of the mean-field Hamiltonian . 
This approach has been initiated by Affleck and Marston in 



Ref. 28 and has been used both in SU(2) and SU(N) mod- 
els. Quite interestingly, in some of the mean-field solutions, 
the product flXij around a plaquette is a complex number 
(X e**^, as if the fermions were picking up a phase due to a 
magnetic field with flux $ threading through the plaquette (in 
convenient units). The finite flux can considerably change the 
band structure of the hopping hamiltonian as well as the cor- 
relations. 

In particular, Hermele et al in Ref. 33 pointed out that 
time-invariance breaking chiral solutions with a uniform flux 
$ are good ground state candidates in a particular large-N 
limit on the square lattice. 

On the honeycomb lattice, the mean-field method has been 
used for the SU(6) symmetric Heisenberg model in Ref. 34, 
where several mean-field solutions have been put forward as 
candidates for the ground state. The lowest mean-field en- 
ergy solution turned out to be a chiral one, with a finite flux 
$ = 27r/3 per hexagon, in line with the ideas put forward 
in Ref. 33. Apart from that, hexamerized solutions with real 
Xi,j hoppings were also found. While the mean-field solu- 
tions give a very useful insight into the possible nature of the 
ground state, they describe free fermions where the on-site 
occupancy is fixed to integer (one if the fundamental represen- 
tation of the SU(N) is considered) on the average only. The 
treatment of the charge fluctuations beyond mean field is quite 
involved. 

There is a complementary approach in which one does not 
search for mean-field solutions, but one takes the ground- 
state wave function of some free-fermion Hamiltonian, and, 
by applying a Gutzwiller projection, one ensures that the oc- 
cupation of each site is exactly one. The projection is done 
numerically, using Monte-Carlo importance sampling, where 
one can efficiently sample the wave function and calculate 
the energy and correlations. ^^'^^ For SU(N) Heisenberg chains 
this approach gives excellent energies and even the correlation 
functions are reproduced with high accuracy.^^"^^ Regarding 
two-dimensional Heisenberg models, the method has been 
applied to the SU(4) model on the square^^ and honeycomb 
lattice,^ ^ and to the SU(3) model on the triangular lattice. As 
for the mean-field solution, introducing a nonzero flux for the 
elementary plaquettes of the hopping Hamiltonian can drasti- 
cally change the band structure, and, after Gutzwiller projec- 
tion, lead to an energy lower than that of a purely 0-flux state. 
This happens for instance in the case of the SU(4) model on 
the honeycomb lattice,^ ^ where introducing a 7r-flux leads to 
an algebraic color liquid that breaks neither the space-group 
nor the SU(4) symmetry. The energy of the 7r-flux state is 
considerably lower than that of the 0-flux state, in which the 
hoppings are all equal. Introducing the 7r-flux leads a band- 
structure with a Dirac-node at the Fermi-level for the quarter 
filled bands (for N colors of the SU(N) model the filling is 



B. Candidate ground states 

Since the family of flux states is infinite, it is of course 
impossible to make a systematic study, and one has to make 
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FIG. 7: (Color online) (a) The 000-flux and (b) the OTTTr-flux Kekule 
states, with 6 sites and 3 hexagons in the unit cell, (c) The tttttt- 
flux and (d) the vrOO-flux states have 12 sites in the primitive unit 
cell (shown by dashed rectangle) or 24 sites in the hexagonal unit 
cell. These states are characterized by two different absolute values 
of hopping amplitudes, td and th. The hopping amplitudes on the 
thin black and grey bonds are td and —td, while on the thick dark 
and light purple bonds the hopping amplitudes are th and —th, re- 
spectively. If the signs of all the hoppings are equal, the product 
of the hoppings around any of the hexagons is positive: this is the 
000-flux state shown in (a). In the other cases, the product of the 
hoppings around some of the plaquettes is negative, corresponding 
to a TT-flux on these plaquettes. (e) The chiral $$$-flux state with 
$ = 27r/3. The red bonds with arrows denote hoppings t/^e^^^^^ 
with complex amplitudes. Here one can also introduce the Kekule 
modulation for the hoppings by changing the sign on the bonds cross- 
ing the boundaries of the unit cell {td bonds), resulting in a 
flux configuration, with ^' = Stt/S. Also and flux 

configurations can be created by introducing complex hoppings to 
the TTTTTT and ttOO case, (f) Brillouin zone of the honeycomb lat- 
tice (black hexagon) with the high symmetry points F = (0,0), 
K = (27r/3^,27r/3), and M = (0,27r/3). The Brillouin-zone 
of the 000, Otttt, and $$^$^-flux states is shown by the dark 

red hexagon, with the high symmetry points Ko = (0,47r/9) and 
Mo = (27r/3A/3, 27r/6). The dark green hexagon stands for the 
Brillouin-zone of the tttttt, ttOO, , and $^$$-flux states with 

= (0,27r/9) and = (7r/3V3, 7r/6). The index in M and 
K refers to the flux of the central hexagon realized by real hopping 
amplitudes. The nearest-neighbor distance is chosen to be unity. 



choices guided by simplicity, intuition, or previous results ob- 
tained on similar models. In the present case, we have decided 
to concentrate on three types of states: 

1. The two Kekule-likc states^^ with flux in a central 
hexagon and flux or tt in the adjacent hexagons called 
respectively 000 and OTTTr-flux states. These states are 
compatible with a 6-site unit cell. In this gauge, the 
hopping amplitudes around the central hexagon are set 
to th, while they alternate between th and td as one goes 
around the two remaining hexagons in the unit cell, as 
shown in Figs. 7(a)-(b). The motivation to study these 
states comes from the on-site and bond color correla- 
tions reported in the tensor network simulations'^^ and 
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in the linear flavor-wave theory,"^^ as weU as our iPEPS 
calculation (Fig. 2) and ED results. 

2. The two Kekule-likc states with tt flux in a central 
hexagon and flux or tt in the adjacent hexagons called 
respectively ttOO and TTTTTr-flux states. The realization 
of these states requires a larger unit cell, as shown in 
Figs. 7(c)-(d). These states are motivated by the results 
of the SU(4) case on the honeycomb lattice.^ ^ 

3. The uniform chiral <l><l><l>-flux state with <l> = 27r/3 
per hexagon (Fig. 7(e)), following the mean-field re- 
sults for the SU(6) Heisenberg model on the honeycomb 
lattice,^^'^^ as well as the uniform ^'^'^'-fiux state, 
where = Stt/S. Both uniform flux states can be 
modulated to achieve a ^^^^^ and a ^^^^ flux config- 
uration. The states with flux <l> in the central hexagon 
can be realized in the 6-site unit cell, while the states 
with flux ^ necessitate a 12 site primitive unit cell (i.e. 
24 site hexagonal unit cell). 

We use the notation t instead of the x ^ distinguish the 
hopping amplitudes set by hand in the variational approach 
from the solutions x of the mean-field approach. Since the 
lattice is bipartite, only the relative sign of the hoppings th 
and td is important, so we can choose > and parametrize 
our results by the ratio td/th- Note that changing the sign of 
the ratio td / allows us to introduce an additional flux tt into 
the hexagons surrounding the central hexagon. Except for the 
chiral phases, both and td can be chosen to be real numbers. 



C. Variational Monte Carlo results 

In the following, we calculate the expectation value of the 
exchange Hamiltonian in the Gutzwiller projected wave func- 
tion using Monte-Carlo importance sampling as a function of 
td for the different realizations of the fluxes. We choose sys- 
tem sizes that have the full symmetry of the hexagonal lat- 
tice and that are compatible both with the 18-site unit cell of 
the SU(3) symmetry broken state with the long-range order 
shown in Fig. 1(a) and with the 24-site unit cell of the tttttt 
and ttOO flux states. This leaves us with two families of clus- 
ters that have 72, 2^ x 72 = 288, 3^ x 72 = 648,... and 
216, 2^ X 216 = 864,... sites. We have used the 72-site clus- 
ter to calculate the energy in a wide parameter range, and we 
have refined the calculation on the 288-site cluster around the 
minima found on the 72-site cluster. In our Monte-Carlo sam- 
pling an elementary update was the exchange of two randomly 
chosen fermions at arbitrary sites with different colors. Each 
run had 10^° (2 x 10^°) elementary updates for the 72 (288) 
size cluster. To avoid any autocorrelation effect, we performed 
measurements after every 1000 (10000) updates, which corre- 
sponds to about 5 times the autocorrelation time. We averaged 
over all the bonds when calculating the energy, and the errors 
of the energy /site values are of the order of 10~^, thus much 
smaller than the symbol sizes on the plots. Furthermore, the 
boundary conditions (periodic or antiperiodic) of the hopping 
hamiltonian were chosen to avoid the degeneracy of the free- 
fermion ground state, if possible at all. 
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FIG. 8: (Color online) (a) Nearest-neighbor bond energy as a func- 
tion of td/th in the 72 site cluster. As we change the sign of td/th 
from negative to positive, we shift between the OTTTr-flux and 000- 
flux states (green squares), or between the ttOO-Aux and TrTrvr-flux 
states (purple circles). The energy of the chiral $$$-flux state with 
$ = 27r/3 is compared to the td = th case, while the energy of 
the chiral $$^$^-flux state is compared to the td = —th case. 
The $^$^$^-flux and $^$$-flux states have a higher energy (not 
shown). The inset of (a) shows the energies around td/th = —1 for 
the 72 and 288 site clusters. The ground state is degenerate for the 
TTTTTT-flux statc whcn td/th > 1, and this is the origin of the scat- 
tered energy values, (b) The energies of the d and h bonds versus td- 
The hexamerization ({Ph) < (Pd)) is more extended for the ttOO- 
flux state than for for the OTTTr-flux state, (c) Schematic drawing of 
the extension of the hexamerized (plaquette) and dimerized phases 
that can be read off from the bond energies given in (b). The arrows 
denote the minima of the energies shown in (a). 



The Monte-Carlo results for the energy per site, given as 
E = (Ph) + (^(i)/2, are shown in Fig. 8(a), while (Ph) 
and (Pd), the expectation values of the exchange operator on 
the bonds of the hexagons and on the dimers, respectively, 
are shown in Fig. 8(b). For td = 0, the wave function is a 
product of decoupled hexagons. For a single hexagon with 
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FIG. 9: (Color online) The band structures of the free fermion Hamiltonian along the path MT^rKTrMj^ in the Brillouin-zone shown in Fig. 7(c) 
for (a) the hexamerized ttOO-Aux state (td = —th), (b)-(d) the dimerized OTrvr-flux state for different values of td- For td > 1 the Fermi surface 
is a Dirac point at Kq — the minimal energy of the corresponds to the case when the Fermi sea touches the s/\th\ = — 1, F point for td = —th 
[plot (c)]. The dark and light purple lines denote occupied and empty bands, respectively, the green dashed line shows the Fermi energy. 



TT-flux (antiperiodic boundary conditions), the variational cal- 
culation leads to (Ph) = -23/30 = -0.7667, very close 
to (Ph) = —0.7676 from exact diagonalization of a six site 
SU(3) Heisenberg chain. Since the ground state of a hexagon 
is a singlet, the exchange between the decoupled hexagons 
is (Pd) = 1/3. This is seen in Fig. 8, where the tttttt and 
ttOO states meet attd = 0. The tt phase in the hexagon with 
nonvanishing hopping amplitudes makes the free fermion lev- 
els twofold degenerate with e = ±V3th and 0, so that we 
have a closed shell condition for two fermions, i.e. the two 
fermion ground state is non-degenerate with energy 
(we need to put 2 fermions of each color to reach 6 fermions 
per hexagon, corresponding to a filling of one fermion per 
site). By contrast, the energy levels of a hexagon with 
flux are e = ±2t^ and e = ±th, the latter levels being 
twofold degenerate. In this case the two-fermion ground state 
is also twofold degenerate, with energy —3th- So, decou- 
pled hexagons prefer to have tt flux. Hexamerization ({Ph) < 
(Pd)) is present for < td/th < I in both the tttttt and 000 
phases, and dimerization takes over for 1 < td/th- In the Otttt 
state the hexamerization persists in the —0.5 ^ td/th < re- 
gion, while for the ttOO case the hexamerization is present in 
a larger window, for —1.5 < td/th < 0. For th ^ the value 
of {Pd) tends to -1. 

Among all states that we have investigated, the TTOO-flux 
configuration provides the lowest energy per site, -0.6912 on 
the 72-site cluster. The minimum occurs around td/th ^ — 1, 
where the projected state shows strong hexamerization (the 
(Ph) is significantly smaller than (Pd), see Fig. 8). The min- 
imum of the OTTTT-flux wave function showing dimerization 
also occurs around td/th ^ —1, with an energy -0.6807 per 
site (TVs = 72) that is slightly higher than that of the hexam- 
erized TTOO-flux state. Furthermore, the energy of the dimer- 
ized state depends only weakly on the cluster size, whereas 
the energy of the hexamerized state is significantly lowered 
when going to the 288-site cluster [cf. inset in Fig. 8(b)]. This 
shows that in the thermodynamic limit the hexamerized TTOO- 
flux state has clearly a lower variational energy than the dimer- 
ized OTTTT-flux state, in agreement with the iPEPS results. The 



energy of the tttttt, 000, and chiral states are all above the ttOO 
and Otttt states.^ 

Next, let us investigate some additional properties of the 
ground state, the hexamerized TTOO-flux state, and of its main 
competitor, the dimerized OTTTT-flux state. It is quite interest- 
ing that the minimum of the energy is around td = —th in both 
cases. We have no explanation why this is so for the hexam- 
erized TTOO-flux state since there is nothing particular happen- 
ing at that point in the free fermion band structure (Fig. 9(a)), 
the occupied bands being well separated from the unoccupied 
bands. By contrast, the Fermi level for the dimerized OTTTT- 
flux state is inside the bands (Fig. 9(b)-(d)), and the td = —th 
point is a special one where the Fermi energy sits both at a 
band edge and at a Dirac point: it separates a fermionic state 
with a finite Fermi surface from a state where the Fermi sur- 
face reduces to a high symmetry point Kq, at which there is a 
Dirac-cone. 

The differences between the bond-energies on the dimers 
and hexamers for the minimal energy TTOO-flux state in the 
VMC is (Pd) - (Ph) ^ 0.6 - 0.7, i.e. larger than the 
iPEPS result of ^ 0.3 - 0.45. For the OTTTT-flux dimerized 
case the VMC result is ^ —0.35(5), while iPEPS provides 
^ -0.26(2). 

The color-color correlations ex (Pij) — 1/3 decay rapidly 
with the distance, as shown in Fig. 10. Since there is a Dirac- 
cone in the free-fermion spectrum of the dimerized OTTTT-flux 
state for td/th < —1, the question whether the correlations 
decay algebraically, and not exponentially as expected for a 
gapped state, is legitimate. Unfortunately, we cannot deter- 
mine unambiguously the nature of the color-color correlations 
in real space, even when using the results from the 648-site 
cluster, the largest we considered. 



We have also compared the energy of different chiral states on the 72 site 
cluster: E(^^'^') = -0.671, E{^^^) = -0.662, E{^'^^) = 
-0.629, and E(^'^'^') = -0.604, where <^ = 27r/3 and = tt + 
^ = 57r/3. 
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FIG. 10: The expectation value of the (P(r)) = {Pij) operator, 
where r is the distance between the sites i and j, for the (a) dimerized 
Otttt and (b) hexamerized TrOO-flux phase, for the 72 and 288 site 
clusters with id — —th- The size of the symbols is proportional to the 
number of bonds having that value of (P(r)). For r ^ oo the value 
of (P(r) tends to 1/3 (denoted by the thin black line), corresponding 
to the expectation value for the exchange between independent spins, 
which shows the absence of the long-range order. 



Next, we look for the signature of the Dirac point in the 
structure factor 



1 ^ V 

Z ( ^^0'^) - ^ ) cosk • (r^- - ro) , 



(7) 



where the summation is over the Ng sites of the cluster, and 
Ti is the position of site i. The pref actor is chosen so that the 
X^IcgbZa *^(^) ~ ^^^^ satisfied, where the sum is 

over the 3Ns/2 k vectors of the Brillouin zone (with the high 
symmetry points and K^) of the underlying triangular 
lattice which, in addition to the sites of the honeycomb lattice, 
also includes the centers of the hexagons. We find that the 
behavior of the structure factor S'(k) is remarkably different 
for the Otttt and ttOO cases [Fig. 1 1] close to the -point: in 
the former one 5'(k) is peaked in the second Brillouin zone, 
while the latter one is smooth. The position of the peak, when 
folded back to the Brillouin-zone of the 6-site unit cell, is at 
Kq, the wave vector that corresponds to the scattering between 
the Dirac points. 



D. Check of the stability towards a color-ordered state 

Following the iPEPS results, which point to a possible 
SU(N) symmetry breaking in the dimerized state, we further 
investigate the possibility of the formation of long-range or- 
der in the Gutz wilier-projected variational approach as well. 
To this end, we allow for color dependent hopping parameters 
and onsite energy terms, starting from the anticipated mini- 
mum energy variational state with td = —th- For each colored 
fermion we introduce a negative onsite energy eLRo according 
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FIG. 1 1 : (Color online) The spin structure factor 5'(k) in the /c-space 
for the (a) hexamerized vrOO-flux and (b) dimerized OvrTr-flux state, 
calculated from the 648-site cluster, with td = —th- The light-green 
hexagons indicate the Brillouin zone of the honeycomb lattice, with 
the high-symmetry points F, K, and M. The 5'(k) for the Otttt- 
flux state is peaked one-third way on the line between the points 
Ma and in the second Brillouin-zone of the honeycomb lattice, 
as can also be seen in (c), where the ^^(k) is shown along the path 
FM^i^AF for the 288 and 648-site cluster. The S{k) for the ttOO- 
flux phase (open symbols), also with td — —th, shows no singular 
behavior. The spin structure factors for the two flux states are nearly 
identical away from the Ma point. 



to the long-range ordered pattern of Fig. 2(a) to enforce the 
SU(N) symmetry breaking in the Gutzwiller projected wave- 
function. Furthermore, we also allow for the modification of 
the hopping amplitude (denoted by ^lro) for the colors that 
dominate the two sites of the bond, while leaving the hop- 
ping amplitude unchanged for the third color. The sign of 
the hopping parameters are set to display the appropriate flux 
state. Note that upon reversing the sign of ^lro the fluxes on 
the hexagons do not change since, for each color, we change 
the sign of two (or none) of the hopping parameters around 
a hexagon. Fig. 12 shows the energy of the Gutzwiller pro- 
jected state as a function of eLRo and ^lro for the Otttt and ttOO 
flux states. It can be clearly seen that the energy is minimal 
for eLRo = and ^lro = ^/i, i-e. for the case where there is 
no long-range order in the system. We have repeated the cal- 
culation also for the 000 and tttttt case as well, starting from 
td = th, and we found that the long-range color ordered phase 
is stabilized for the 000-flux state with energy E = —0.610 
per site (for tLRo/|^if| = 1-3 and eLRo/|^/i| = -0.4), much 
higher than the lowest energy ttOO solution. 
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FIG. 12: (Color online) The stability of the dimerized and hexamer- 
ized states versus the formation of long range order, in the (a) Otttt 
and (b) ttOO flux state, respectively. For both cases the energy is 
minimal for clro = and tLRo = th, where there is no long range 
order. The calculations were made with stepsizes AtLRo = 0.1|t^| 
and AeLRo = 0.2|t/i|. For tLRO = the fermionic band structure 
collapsed to a few highly degenerate bands, this caused difficulties in 
the Monte Carlo code we used for the calculation. 



V. DISCUSSION AND CONCLUSION 

In this work we showed that the ground state of the SU(3) 
Heisenberg model on the honeycomb lattice has plaquette or- 
der and does not break SU(3) rotational symmetry in color 
space, in agreement with the previous tensor network study in 
Ref. 42. Extrapolations of the plaquette order parameter to the 
infinite D limit reveals that the ground state indeed exhibits 
true long-range order. 

This result is in conflict with the prediction by LFWT from 
Ref. 41 . However, by performing a systematic study of the so- 
lution as a function of the bond dimension D in iPEPS we can 
understand the LFWT prediction as a low-entanglement solu- 
tion which is energetically favorable for small D, but which 
is unstable upon taking more quantum fluctuations into ac- 
count by going to large D. This situation is reminiscent of the 
SU(4) Heisenberg model on the square lattice, where LFWT 
and iPEPS with D = 2 predict a plaquette-color ordered state, 
however, for D > 5 a. dimer-Neel ordered state is stabilized.^^ 



Thus, iPEPS is an ideal tool to compare competing states. 
Unlike standard variational methods, one can study how the 
individual energies of the competing states change upon in- 
creasing the amount of quantum fluctuations in a systematic 
way, by varying the bond dimension. Such systematic anal- 
ysis will be important also for future tensor network studies, 
e.g. for the t-J model where a uniform d-wave superconduct- 
ing state is competing with a striped state at very low ener- 
gies.49'60 

Next, comparing the energy of Gutz wilier projected wave 
functions, we have found that the competition between pla- 
quette formation (hexamerization) and dimerization is also 
present in the variational treatment. We have identified that 
the two competing states originate from a free-fermion wave 
function with ttOO and Otttt fluxes on the hexagons in the 
unit cell, with very different nature: the lower energy ttOO- 
flux state describes a gapped, plaquette (hexamerized) state, 
while the higher energy Otttt flux state is gapless and dimer- 
ized, with a free-fermion wave function having a Fermi-point 
at the Dirac cone. This difference is exemplified in the struc- 
ture factor, in the former case it is a smooth function, in the 
latter it has a peak related to the position of the Dirac cones. 
A dimerized solution with a broken SU(3) symmetry has also 
been identified in the variational treatment which, however, 
has a much larger energy than the two competing states. 

Finally we note that the plaquette ground state is very dif- 
ferent from the one obtained for the same model on the square 
lattice, which exhibits three- sublattice Neel order. ^^'^^ In prin- 
ciple, length-6 loops could also be formed on the square lat- 
tice. However, the energy cost to introduce another weak bond 
(across the hexagons which results in a square lattice) is too 
high, which makes the three- sublattice Neel ordered state en- 
ergetically favorable in this case. 
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